Supertransient magnetohydrodynamic turbulence in Keplerian shear flows 
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A subcritical transition to turbulence in magnetized Keplerian shear flows is investigated using a 
statistical approach. Three-dimensional numerical simulations of the shearing box equations with 
zero net magnetic flux are employed to determine the transition from decaying to sustained turbu- 
lence as a function of the magnetic Reynolds number Rm. The results reveal no clear transition to 
sustained turbulence as the average lifetime of the transients grows as an exponential function of 
Rm, in accordance with a type-II supertransient law. 
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There has been much progress in recent years on the 
understanding of the nature of subcritical transitions to 
turbulence in hydrodynamic flows, wherein the laminar 
state is linearly stable for all Reynolds numbers and fi- 
nite amplitude perturbations are necessary to drive the 
system towards a turbulent state. Numerical and ex- 
perimental works indicate that turbulence in systems 
like plane-Couette and pipe flows is a transient phe- 
nomenon [TJ [5], although this is still a disputed topic 
[31 11] . Although much has been accomplished for hydro- 
dynamic turbulence in this area, the same cannot be said 
about magnetohydrodynamic (MHD) turbulence in sys- 
tems with subcriticality, for which statistical studies are 
scarce. 

A classical problem of subcritical transition to turbu- 
lence is provided by accretion disks, where, in the absence 
of magnetic fields and when shear is the only stratifica- 
tion, Keplerian flows are linearly stable 5 . The rate of 
mass accretion in the disk depends on the outward trans- 
port of angular momentum and one would think that tur- 
bulence is the natural candidate for transport. However, 
both numerical [BHS] and laboratory '5] experiments have 
revealed that, if a subcritical hydrodynamic transition to 
turbulence was to exist in Keplerian flows, then the asso- 
ciated turbulent transport of angular momentum is very 
weak and unable to match observational constrains. 

When a mean vertical fleld is applied, the magnetoro- 
tational instability [TUj provides a linear instability that 
may result in a strong MHD turbulence capable of gen- 
erating adequate rates of angular momentum transport 
[m [12]. If there is no imposed background magnetic 
field, the original (seed) field can decay with time due 
to finite resistivity. In that case, some type of nonlinear 
dynamo process is necessary to amplify and sustain the 
seed field and, consequently, the instability in zero net 
flux simulations is nonlinear |13j . Rincon et al. |14j iden- 
tifled a self-sustaining dynamo solution in a magnetized 
rotating plane-Couette flow in the Keplerian regime, but 



the steady solutions were found only in a narrow range 
of low values of the Reynolds number. Other works have 
reported the existence of nonlinear dynamos in turbulent 
MHD regimes [13], yet, little is known about how turbu- 
lence evolves as a function of the magnetic Prandtl num- 
ber Pm (the ratio of viscosity to resistivity) and whether 
it can be sustained for long runs or not. Recently, Fro- 
mang et al. |15j studied the onset of turbulence in accre- 
tion disks using an unstratified shearing box with zero 
net magnetic flux and found that turbulence disappears 
when the magnetic Prandtl number falls below a criti- 
cal value that seems to be a decreasing function of the 
kinetic Reynolds number. Note, however, that vertical 
stratification might change this picture |16j . 

The present letter explores the onset of turbulence in 
three-dimensional, zero net flux, MHD simulations of Ke- 
plerian shear flows as a function of the magnetic Reynolds 
number Rm. We focus on the statistical aspects of the 
subcritical transition, which is an approach that has been 
unexplored in MHD simulations, despite its popularity in 
the hydrodynamics community. From our results, there 
seems to be no clear transition to sustained turbulence. 
Instead, the average lifetime of the turbulence r grows 
as an exponential function of Rm, following a type-II 
supertransient law |17j . 

We solve the shearing box equations |18| , which repre- 
sent the dynamics in a box moving with the disk's angular 
velocity ^(tq), where rg is a flducial radius that deter- 
mines the location of the center of the box. In Cartesian 
coordinates, with (/)—>■ y and r — > tq + x, the shearing 
box equations read as 



(9fV + V • Vv = — VP H (V X B) X B 
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V • B = 0, V • V = 0, (3) 

where V, — r^^^^ for a Keplerian disk, S = —rdr^ = 
(3/2)r2, r] is the magnetic difFusivity, fiQ is the magnetic 
permeabihty, p is the gas density, and P is the pressure. 
The fluid velocity can be decomposed as v = vq + u, 
where the steady state solution is given by the shear 
flow Vq = —Sxy and u is the perturbation field. The 
equations for u are solved with the pseudospectral code 
described by Lesur et al. [3 [19] using shearing sheet 
boundary conditions [TT] . The xyz aspect ratio employed 
is 1 X TT X 1. The kinetic and magnetic Reynolds numbers 
are defined as Re = S(P /v and Rm = ScP/r], respec- 
tively, where d is the shearwise size of the box. Time is 
measured in shear time units, 1/S. 

Transition to turbulence is investigated by fixing Re 
while varying Rm. The value of Re is taken from Fig. 
11 of Fromang et al. [TS], where for Re = 3125 there is a 
transition to turbulence somewhere between Rm = 6250 
and Rm — 12500, thus, the magnetic Prandtl number 
Pm = Rm/Re > 1. Since we are not interested in the 
kinematic dynamo regime, the initial growth of a small 
seed magnetic field is not important to our discussions. 
Our question is: given a turbulent state, can it be sus- 
tained for all time? If not, how long it takes, on aver- 
age, to decay? The turbulent states are obtained from 
long runs where the velocity and magnetic fields are ini- 
tialized with a large-scale random noise (only the four 
largest modes are nonzero at t = 0) with amplitude of 
the order of 1. After some time, energy spreads to other 
modes and the random fields settle to a turbulent state. 
Figure [T] illustrates the typical turbulent patterns found 
at Re = 3125 and Rm = 11000 for the radial component 
of the velocity (upper panel) and magnetic (lower panel) 
fields at t = 500. Since Pm > 1 and there is no proper 
large scale dynamo, the characteristic spatial scales in 
the magnetic fields are smaller than in the velocity fields. 
The numerical resolution employed is 64 x 128 x 64, which 
is the same resolution used for this problem in Ref. [TS]. 
The convergence of the code is checked by comparing the 
time-averaged power spectra obtained with this run with 
a higher-resolution run (128 x 256 x 128). The results 
are shown in Fig. [2j where the solid black line repre- 
sents the lower-resolution kinetic energy spectrum and 
the solid red line the lower-resolution magnetic energy 
spectrum. The higher-resolution spectra are plotted in 
dashed lines and agree quite well with the lower resolu- 
tion runs. There is a disagreement in the small scales 
in the kinetic spectra, but this has a low impact on the 
global behaviour. Table 1 contains time-averaged val- 
ues for the kinetic and magnetic energies, as well as for 
the Reynolds and Maxwell stresses, which are defined, 
respectively, as a^ey = (uxUy) and aMax = -{BxBy). 
The quantities are computed with both resolutions for 
Re = 12500. In the remainder of this letter, we adopt 




FIG. 1. (Color online) Volume rendering plots of Ux (a) and 
Bx (b) for a turbulent state at Re = 3125 and Rm = 11000. 

the lower resolution. Note that the lack of a properly 
defined inertial range in the power spectra indicates that 
the turbulence is not fully developed. 

TABLE I. Time-averaged quantities for Rm — 12500. 



Res. 


64 x 128 x 64 


128 X 256 X 128 


{Ekin)i 
{o!Max)f- 


2.11 X 10"^ 
8.45 X 10~^ 
5.22 X 10"* 
3.68 X 10~^ 


2.07 X 10"^ 
8.52 X 10"^ 
5.24 X 10"* 
3.69 X 10"-^ 



The turbulence shown in Fig. [T] eventually decays to- 
ward the laminar state. The decay time is a function of 
the initial perturbation and Rm. Figure [3] shows an ex- 
ample of a time series of the magnetic (red) and kinetic 
(black) energies for Re = 3125 and Rm = 12500. Af- 
ter more than 6000 time units the turbulence suddenly 
decays and the system converges to the laminar state. 

The high variability of the turbulence decay times as 
a function of the initial condition for a given Rm sug- 
gests the presence of a chaotic saddle in the phase space 
[H 1^ . The role of chaotic saddles in generating chaotic 
transients in dynamical systems is well known and, as 
mentioned in previous studies on spatiotemporal systems 



3 




Or 
-1 

3 -2- 

M .3- 
O 



-T- ▼ 



• 


Rm = 


9000 


■ 


Rm = 


9375 


♦ 


Rm = 


9750 




Rm = 


10000 


▼ 


Rm = 


10500 


* 


Rm = 


11000 



200 400 600 800 1000 
t 



FIG. 2. (Color online) Time-averaged spectra of kinetic 
(solid black line) and magnetic (solid red line) energies for 
Re = 3125 and Rm = 11000. The dashed lines represent the 
spectra computed with doubled resolution. The k~^^^ dotted 
line is a guide to the eye. 



FIG. 4. (Color onUne) The logarithm of P{t), the probabihty 
of turbulent lifetime > i, as a function of t for Re = 3125 and 
six different values of Rm; the straight lines are fitted with 
P(t, Rm) = exp[— t/r(_Rm)], where r is the characteristic tur- 
bulence lifetime. 
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FIG. 3. (Color online) Time series of the magnetic (red) 
and kinetic (black) energies showing transient turbulence for 
Re = 3125 and Rm = 12500. 



|21j . the average lifetime of the chaotic transients grows 
as a function of the fractal dimension of the stable man- 
ifold of the chaotic saddle. In order to obtain the aver- 
age lifetime r(i?m) of the turbulent transients we com- 
pute P{t), the probability of finding a turbulent state at 
time t. For each value of Rm, Pit) is obtained from a 
set of more than 50 and up to 90 different simulations. 
Each simulation has a turbulent state as initial condi- 
tion and the simulation stops when either t = 1000 or 
the laminar state is reached, generating a data set of 
relaminarization times. The turbulent initial conditions 
are taken from a long run at Rm = 13000. To ensure 
that the initial conditions represent uncorrelated states, 
they were selected from the referred long turbulent run 
with 500 shear-time units apart from each other. Fig- 
ure [4] shows the logarithm of P(t) as a function of t for 
Re = 3125 and a range of different values of Rm. The 
straight lines are least square fits that represent the re- 
lation P{t,Rm) = exp[— i/r(i?TO)], which is expected for 
transients due to chaotic saddles |21 [12] • The values of 
1/r are readily obtained from the slopes of the fitted 
lines. As Rm increases, the slopes decrease, refiecting 



longer decay times. The convergence is not very good for 
the lowest values of Rm, probably due to a lower fractal 
dimension of the stable manifold of the chaotic saddle, 
which implies that a much higher number of initial con- 
ditions would be necessary to obtain better statistics. 

A graph of r versus Rm is plotted in Fig. [5j The top 
panel shows the data in linear scales and reveals that the 
inverse lifetime can be well fitted with an exponential law 
as a function of Rm. The dashed curve is given hy 1/t — 
exp(5.3 — 0.00118i?m). The lower panel shows the same 
data in log-lin scales. These results indicate that the 
turbulence lifetime follows a type-II supertransient law 
|17j as a function of Rm. In this case, there is no apparent 
transition to sustained turbulence, and the laminar state 
is the global attractor. If there were a critical Rm for 
transition to sustained turbulence via a crisis mechanism, 
a linear decay of 1 /r with Rm would have been expected 
[3]. The results for a linear fit using only the last four 
values of Rm are given by the dotted blue line in Fig. 
[sja), which represents 1/r 0.0144 - 1.273 x IQ-^ Rm. 
This fit provides a critical point Rm^ ^ 11312, which is 
clearly incorrect, since Fig. [3] shows that even at Rm = 
12500 the turbulence is still decaying. 

Note that the high values of Re and Rm shown in this 
letter are based on the shear rate S. If the Reynolds num- 
bers are computed based on the root mean square veloc- 
ity, one obtains a maximum of Re ~ 180 and Rm ~ 650. 
Thus, it has been shown that, at least for intermedi- 
ate Reynolds numbers and magnetic Prandtl numbers 
greater than one, MHD turbulence in the incompressible 
shearing box equations with zero net magnetic fiux is not 
self-sustained. The average decay time r follows a type-II 
supertransient law, with an exponential increase of r as a 
function of the magnetic Reynolds number. This is in ac- 
cordance with results found for subcritical transitions to 
turbulence in hydrodynamic fiows such as the pipe flow 
[5] and plane Couette flow [1 , where turbulence seems to 
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FIG. 5. (Color online) The characteristic turbulent lifetimes 
as a function of Rm, obtained from fits to the data shown in 
Fig. [4] (a) linear scales, the dashed curve is an exponential 
fit and the dotted straight line is a linear fit using only the 
last four points; (b) same data as in (a) but in log-lin scales, 
with the exponential fit. 



be always transient. Supposing that our conclusions can 
be extended to realistic values of the Reynolds numbers, 
this suggests that the transient dynamics is much more 
relevant to the accretion problem than the final attractor. 
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